using System;
using System.Linq;

using Atomic.Structures;
using Atomic.Thermodynamics.StateEquations;

namespace Atomic.Thermodynamics
{
	public static class ThermodynamicsExtensions
	{
		public static double PressureNumerical(this IStateEquation stateEquation, double v)
		{
			double dv = 0.0001;
			double f_v = (stateEquation.Energy(v + dv) - stateEquation.Energy(v)) / dv;

			return -f_v;
		}

		public static double BulkModulusNumerical(this IStateEquation stateEquation, double v)
		{
			double dv = 0.0001;
			double p_v = (stateEquation.Pressure(v + dv) - stateEquation.Pressure(v)) / dv;

			return -v * p_v;
		}

		public static double PressureNumerical(this IHelmholtzPotential potential, double t, double v)
		{
			double dv = 0.0001;
			double a_v = (potential.FreeEnergy(t, v + dv) - potential.FreeEnergy(t, v)) / dv;

			return -a_v;
		}

		public static double PressureTemperatureDerivativeNumerical(this IHelmholtzPotential potential, double t, double v)
		{
			double dt = 0.0001;
			double p_t = (potential.Pressure(t + dt, v) - potential.Pressure(t, v)) / dt;
			double a_vt = -p_t;

			return -a_vt;
		}

		public static double BulkModulusNumerical(this IHelmholtzPotential potential, double t, double v)
		{
			double dv = 0.0001;
			double p_v = (potential.Pressure(t, v + dv) - potential.Pressure(t, v)) / dv;
			double a_vv = -p_v;

			return v * a_vv;
		}

		public static double PressureNumerical(this IHelmholtzGrandPotential potential, double t, double v, Composition composition)
		{
			double dv = 0.0001;
			double a_v = (potential.FreeEnergy(t, v + dv, composition) - potential.FreeEnergy(t, v, composition)) / dv;

			return -a_v;
		}
		
		public static ChemicalPotential ChemicalPotentialNumerical(this IHelmholtzGrandPotential potential, double t, double v, Composition composition)
		{
			Atom atom;
			double x;
			UnaryComposition(composition, out atom, out x);

			double dx = 0.0001;

			// Derivative of Helmholtz free energy.
			double a_x = (potential.FreeEnergy(t, v, new Composition(new Constituent(atom, x + dx))) - potential.FreeEnergy(t, v, new Composition(new Constituent(atom, x)))) / dx;

			return new ChemicalPotential(atom, a_x);
		}

		public static ChemicalPotential ChemicalPotentialNumerical(this IGibbsGrandPotential potential, double t, double p, Composition composition)
		{
			Atom atom;
			double x;
			UnaryComposition(composition, out atom, out x);

			double dx = 0.0001;

			// Derivative of Gibbs free energy.
			double g_x = (potential.FreeEnergy(t, p, new Composition(new Constituent(atom, x + dx))) - potential.FreeEnergy(t, p, new Composition(new Constituent(atom, x)))) / dx;

			return new ChemicalPotential(atom, g_x);
		}

		private static void UnaryComposition(Composition composition, out Atom atom, out double count)
		{
			Atom[] atoms = composition.Atoms.Distinct().Where(a => a != Atom.Vacancy).ToArray();

			if (atoms.Length != 1)
			{
				throw new ArgumentException();
			}

			atom = atoms[0];
			count = composition.Count(atom);
		}
	}
}
